#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "../common/general.h"
#include "../common/minilog.h"

#define VERBOSE		INFO | RESULT | WARNING | ERROR | CRITICAL
#define LOGGING		INFO | RESULT | PROCESS
#define COMMENT		"This is the canonical version of HS, using Gray-Encoding.\nRandom Generator: Mersenne twister"
#define LOGFILENAME	"HSGC.hsgc"
#define ALGORITHM	"HSGC"

/********************************************************************************
	Harmony Search Parameters
********************/
#define max_iterations		500000
#define hms					10
#define hmcr				0.95
#define par					0.005
#define fw					4
#define test_function		7
#define precision_digits	8
#define sample_step			1e-10
#define eval_stamp			100000

/*
#define dimension			100
#define lower_bound			-1000.0
#define upper_bound			+1000.0
*/
#include "../benchmark/group02fn.h"

/********************************************************************************
	Harmony Search Global Variables
********************/
struct solution{
	char* bit_stream;
	double fitness;
};
struct solution* hm;
int nvarbits;
int nvecbits;
double* power_two;
double* last_x;
double precision;
int eval_counter;
char unique_logfilename[256];
#include "../common/common.h"
#include "../common/encoding.h"


/********************************************************************************
	Subordinate Functions
********************/
int solComp(const void* a,const void* b) {
	double pa,pb;
	pa=((struct solution*)a)->fitness;
	pb=((struct solution*)b)->fitness;
	if (pa>pb) return 1;
	else if (pa<pb) return -1;
	return 0;
}
void copy(struct solution from,struct solution *to){
	int i;
	to->fitness=from.fitness;
	for(i=0;i<nvecbits;i++)
		to->bit_stream[i]=from.bit_stream[i];
}

void copy_var(struct solution from,struct solution *to,int svar,int dvar){
	int i,db,sb;
	db=dvar*nvarbits;
	sb=svar*nvarbits;
	for(i=0;i<nvarbits;i++)
		to->bit_stream[db+i]=from.bit_stream[sb+i];
}

/********************************************************************************
	Harmony Search Main Functions
********************/
void improvise(int index){
	int i;
	for (i=0;i<nvecbits;i++)
		hm[index].bit_stream[i]=randint3(2);
}
void improvise_var(int index,int var){
	int i,b;
	b=var*nvarbits;
	for (i=0;i<nvarbits;i++)
		hm[index].bit_stream[b+i]=randint3(2);
}

void eval(struct solution *sol){
	vector_gray2dec_range(sol->bit_stream,last_x);
	sol->fitness=func(last_x);
	eval_counter++;
}

int init(){
	int i;
	eval_counter=0;
	build_min_vector();
	if(precision_digits)
		precision=1.0/pow(10,precision_digits);
	else precision=0.0;
	cal_nvarbits();
	nvecbits=dimension*nvarbits;
	hm=(struct solution*) malloc(sizeof(struct solution)*hms*2);
	power_two=(double*) malloc(sizeof(double)*(nvarbits+1));
	last_x=(double*) malloc(sizeof(double)*dimension);
	power_two[0]=1.0;
	for(i=1;i<=nvarbits;i++) {
		power_two[i]=power_two[i-1]*2.0;
	}
	for(i=0;i<hms*2;i++){
		hm[i].bit_stream=(char*)malloc(nvecbits);
		improvise(i);
		eval(&hm[i]);
	}
	qsort(hm,hms*2,sizeof(struct solution),solComp);

}
void harmony_consideration(){
	int i,j;
	for(i=0;i<hms;i++){
		for(j=0;j<dimension;j++){
			if(random3()<hmcr)
				copy_var(hm[(int)(random3()*hms)],&hm[hms+i],j,j);
			else
				improvise_var(hms+i,j);
		}
	}
}
void pitch_adjustment(){
	int i,j,k,ind;
	for(i=hms;i<hms*2;i++){
		for(j=0;j<dimension;j++){
			for(k=0;k<fw;k++){
				if (random3()>=par) continue;
				ind=randint3(nvarbits);
				if(hm[i].bit_stream[(j*nvarbits)+ind]==0)
					hm[i].bit_stream[(j*nvarbits)+ind]=1;
				else
					hm[i].bit_stream[(j*nvarbits)+ind]=0;
			}
		}
		eval(&hm[i]);
	}
}
int stopping_criterion(){
	double r;
	r=euclid_dist(last_x,min_vector);
	if(r<precision) return 1;
	return 0;
}

/********************************************************************************
	Display and Debug Functions
********************/
void show_info(FILE* out){
	fprintf(out,"/*\n%s\n*/\n",COMMENT);
	L(INFO)fprintf(out,"Logging Filename                 : %s\n",unique_logfilename);
	fprintf(out,"Algorithm                        : %s\n",ALGORITHM);
	fprintf(out,"Test Function                    : %s\n",func_name);
	fprintf(out,"Number of Variables              : %d\n",dimension);
	fprintf(out,"Number of Bits per Variable      : %d\n",nvarbits);
	fprintf(out,"Number of Bits per Vector        : %d\n",nvecbits);
	fprintf(out,"Variables Lower Bound            : %.15f\n",lower_bound);
	fprintf(out,"Variables Upper Bound            : %.15f\n",upper_bound);
	fprintf(out,"Variables Precision              : %.15f\n",precision);
	fprintf(out,"Maximum Iterations               : %d\n",max_iterations);
	fprintf(out,"Harmony Memory Size              : %d\n",hms);
	fprintf(out,"Harmony Memory Consideration Rate: %.15f\n",hmcr);
	fprintf(out,"Pitch Adjusting Rate             : %.15f\n",par);
	fprintf(out,"Fret Width                       : %d\n",fw);
	fprintf(out,"-------------------------------------------------------------\n");
}
void show_memory(){
	int i,j;
	for(i=0;i<hms*2;i++){
		for(j=0;j<nvecbits;j++)
			printf("%d",hm[i].bit_stream[j]);
		printf(": %.7f\n",hm[i].fitness);
	}
}

int main(){
	int n=0;
	int i;
	double last_best_fitness=1e200;
	double time_stamp;
	double dist;
	clock_t start;
/*	srand(time(NULL));*/
	mt_init();
	get_unique_filename(LOGFILENAME,unique_logfilename);
	minilog_start(unique_logfilename,VERBOSE,LOGGING);
	init();
	V(INFO)show_info(stdout);
	L(INFO)show_info(LOGFILE);
	start=clock();
	while(n<max_iterations){
		harmony_consideration();
		pitch_adjustment();
		qsort(hm,hms*2,sizeof(struct solution),solComp);
		eval(&hm[0]);eval_counter--;
		dist=euclid_dist(last_x,min_vector);
//		VP(PROCESS,"[%06d] f=%.15f\td=%.15f             \r",n,hm[0].fitness,dist);
		VP(PROCESS,"[%07d] f=%.15f\td=%.15f          \r",eval_counter,hm[0].fitness,dist);
		if(fabs(hm[0].fitness-last_best_fitness)>=sample_step || eval_counter%eval_stamp==0) {
//			LP(PROCESS,"%06d\t%.15f\t%.15f\n",n,hm[0].fitness,dist);
			LP(PROCESS,"%07d\t%.15f\t%.15f\n",eval_counter,hm[0].fitness,dist);
			last_best_fitness=hm[0].fitness;
		}
		n++;
	}
	time_stamp=(double) (clock()-start)/CLOCKS_PER_SEC;
//	LP(PROCESS,"%06d\t%.15f\t%.15f\n",--n,hm[0].fitness,dist);
	LP(PROCESS,"%07d\t%.15f\t%.15f\n",eval_counter,hm[0].fitness,dist);
	VP(PROCESS,"\n");
	LOG(RESULT,"-----------------\n");

	LOG(RESULT,"Best Fitness           : %.15f\n",hm[0].fitness);
	LOG(RESULT,"Closest Distance       : %.15f\n",dist);
	LOG(RESULT,"Time                   : %05.2fs\n",time_stamp);
	LOG(RESULT,"-----------------\n");
	for(i=0;i<dimension;i++){
		LOG(RESULT,"x[%d]= %.15f\n",i,last_x[i]);
	}
	minilog_stop();
	return 0;
}
